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ABSTRACT 

A new cosmic shear analysis pipeline SUNGLASS (simulated UNiverses for Gravitational 
Lensing Analysis and shear Surveys) is introduced. SUNGLASS is a pipeline that rapidly 
generates simulated universes for weak lensing and cosmic shear analysis. The pipeline 
forms suites of cosmological N-body simulations and performs tomographic cosmic 
shear analysis using line-of-sight integration through these simulations while saving 
the particle lightcone information. Galaxy shear and convergence catalogues with re- 
alistic 3D galaxy redshift distributions are produced for the purposes of testing weak 
lensing analysis techniques and generating covariance matrices for data analysis and 
cosmological parameter estimation. We present a suite of fast medium resolution sim- 
ulations with shear and convergence maps for a generic 100 square degree survey out 
to a redshift of z = 1.5, with angular power spectra agreeing with the theory to better 
than a few percent accuracy up to £ = 10'^ for all source redshifts up to z = 1.5 and 
wavenumbers up to £ = 2000 for the source redshifts z ^ 1.1. At higher wavenum- 
bers, there is a fail ure of the th e oretic al lensing power spectrum reflecting the known 
discrepancy of the ISmith et al.l (|2003[ ) fitting formula at high physical wavenumbers. 
A two-parameter Gaussian likelihood analysis of erg and flm is also performed on the 
suite of simulations, demonstrating that the cosmological parameters are recovered 
from the simulations and the covariance matrices are stable for data analysis. We find 
no significant bias in the parameter estimation at the level of ^ 0.02. The SUNGLASS 
pipeline should be an invaluable tool in weak lensing analysis. 

Key words: Gravitational lensing - Cosmology: large scale structure of Universe - 
Methods: iV-Body simulations 



1 INTRODUCTION 

Cosmic shear analysis is an excellent met hod for prob 



ing the dark Universe (for r eviews, s ee Melliei _ 19991: 



BM^hnani^^_SchMide3^ 



200ll: iRefregieJ l2003l : ISchneideij 



20061 : iMunshi et alJ bOOa : iMassey et aLr[2010r and refer- 
ences therein). It is also a reasonably new field of research 
with cosmic shear first being observed just ten years ago 
llBacon et al.ll2000l: iKaiser et all [2OO0I : I Van Waerbeke etal] 
120001 : IWittman et al. I I2OOOI ). Weak gravitational lensing ef- 
fects on a cosmic scale are a mere 1% change in shape and 
observational systematics makes the measurement of these 
changes challenging. However, the combination of the well- 
understood underlying physics and the expected precision 
of cosmological parameter estimation make the effort worth- 
while. 

Next generation telescope surveys will observe more of 
the sky than ever before and the volume of data they will 
produce is unprecedented. Future surveys promise to deter- 
mine the equation of state of dark energy to 1% as well as 

* E-mail: aak@roe.ac.uk 



probing the possibilities of extra dimensional gravity models 
and alternative cosmologies. The first Pan-STARR^| tele- 
scope is currently undertaking a cosmic shear survey of 
the entire visible sky from its location in Hawaii and new 
projects such as VST-KID^|, DE^, HALCQ, Euchcfl and 
LSSI0 and are planned to perform wide field cosmic shear 
surveys, measuring both large, linear scales, and small, non- 
linear scales. 

Due to the relative youth of this field, techniques are 
still being developed to exploit the weak lensing data from 
these surveys to provide further understanding on the na- 
ture of the Universe. To realise the potential of these new 
telescope surveys and to test new weak lensing analysis tech- 
niques, challenges must be met. To achieve the small sta- 
tistical errors required, experiments require full end-to-end 



^ Pan-STAR RS|http:/ /pan-starrs. ifa.hawaii.edu/"] 

2 VST-KIDS 'http://www.astro-wise.org/ projects/KIDS/ 
^ DES https: / /www. darkenergysurvey.org/ 
* Rhodes et al. , in preparation 
^ Euclid http://sci.esa.int/euc lid] 
LSST ht tp777www.lsst.org/| 
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simulations of huge volumes which also probe the non-linear 
regime to assist in understanding the limitations of the anal- 
ysis techniques. Simulations offer data sets with known pa- 
rameters which are essential when testing analysis pipelines. 
Simulations can also include effects which may be difficult 
to model theoretically, such as source clustering and galaxy 
alignments, as well as other systematics and real- world ef- 
fects. An additional role for simulations is in accurate es- 
timation of the covariance of observable quantities. This 
is needed for the analysis of surveys and analytic approx- 
imati ons can be wholly inadequate (e.g. ISemboloni et al.l 
[2003). Monte Carlo analysis can be performed with simu- 
lations to provide covariance matrices that are required for 
data analysis and cosmological parameter estimation. Sim- 
ulations are also required for rigorous testing and develop- 
ment so all analysis methods can be analysed blindly before 
the same techniques are applied to real data. To address 
these challenges, the SUNGLASS, Simulated UNiverses for 
Gravitational Lensing Analysis and Shear Surveys, pipeline 
has been developed to produce simulations and mock shear 
and convergence catalogues rapidly for weak lensing and cos- 
mic shear analysis. The purpose of this paper is to introduce 
SUNGLASS and show rigorous testing of its outputs. 

Many weak lensing studies use si mulations with very 
high resolution to run their analysis fe.g.lFosalba et al]l2008l: 
Hilbert et all l2009l: iTevssier et al] |2009| : ISchrabback et"aLl 



2OI0I '). The computational cost of running these simula- 
tions is high and consequently there is often only a single 
realisation available. However, it is very important to en- 
sure that covariance matrices calculated from these simu- 
lation s are not contamina ted by correlations in the simula- 
tions (|HartIap et al.|[2007l ). In order to ensure uncorrelated 
data, a Monte Carlo suite of simu lations should be used to 
determine the covariance matrix (|Sato et al.ll2009h . In this 
work, 100 independent simulations were constructed using 
SUNGLASS. 

To date, there are still reasonably few weak lensing sim- 
ulations available. Of the few that are available, many im- 
plement a ray-tracing technique where light rays are prop- 
agated from an observer to a lensing source plane (e.g . 
Jain et all l200d: IVale fc White! l20o3: lHilbert etHI l2009l: 



Tevssier et all I2OO9I: ISato et al.l 120091 : [Dietrich fc Hartlanl 
2010l : IVafaei et al.ll2010l )~ Ray-tracing is computationally in- 



tensive and time consuming when solving the full ray-tracing 
equations. If the Born approximation is used in the ray- 
tracing, the time to run the analysis is reduced but the 
process is still computationally intensive and the simulation 
data still needs to be binned in three dimensions to per- 
form the calculations. An alternative to ray-tracing is line- 
of-sight integration, which uses the Born approximation to 
calc ulate rapidly the we a k lensing signal thro ugh a lightcone 
(e.g. IWhite fc Hull2000l : [Posalba et ahlbOOSl ). This method 
is not suitable in the strong lensing regime but in the weak 
lensing regime, it is rapid and requires fewer computational 
resources than ray-tracing techniques. In this paper, a new 
line-of-sight integration technique, implemented in the SUN- 
GLASS pipeline, for measuring convergences in an N-body 
simulation is introduced. This new method is rapid and can 
be run on a single processor of a desktop computer. In con- 
trast to ray-tracing, the method does not bin in the radial 
direction, using all of the redshift information available. Al- 
though the catalogues are suitable for real-space analysis. 



SUNGLASS analyses and tests our mock weak lensing sur- 
veys in Fourier space, using power spectra, as it is possible 
to cleanly distinguish between linear and nonlinear regimes 
in Fourier space. We are also able to easily identify scales 
where the simulations are reliable by determining the region 
of the power spectrum in Fourier space that lies between the 
size of the simulated volume at low wavenumbers and shot- 
noise due to particle discreteness and pixelization effects at 
high wavenumbers. 

The outline of this paper is as follows. Section [2] intro- 
duces the SUNGLASS pipeline. Details of the simulations are 
in section [27l] and the line-of-sight integration method for de- 
termining shear and convergence without radial binning is 
described in section Section [2.31 presents the shear and 
convergence power spectrum analysis and section [T4l deals 
with the generation of the mock galaxy shear catalogues. An 
application of the mock catalogues is discussed in section [3] 
where Gaussian likelihood estimates of Q.m and erg are per- 
formed. A summary of the pipeline and methods concludes 
the paper in section |3] 



2 DETAILS OF THE SUNGLASS PIPELINE 

SUNGLASS is a pipeline that generates cosmic shear and 
convergence catalogues using N-body simulations. The 
pipeline creates mock galaxy shear catalogues that can be 
used to test the cosmic shear analysis software used on tele- 
scope survey data sets. The nature of the pipeline also al- 
lows many simulation realisations to be generated rapidly to 
produce covariance matrices for data analysis and cosmolog- 
ical parameter estimation. The pipeline begins by creating 
a suite of cosmological N-body simulations. Lightcones are 
generated through the simulations and tomographic shear 
and convergence maps are determined using line-of-sight in- 
tegrations at multiple lensing source redshifts. Finally, mock 
galaxy catalogues with fully 3D shear and convergence in- 
formation and galaxy redshift distributions are assembled 
from the lightcones and the tomographic shear and conver- 
gence planes. The following sections detail each step of the 
SUNGLASS pipeline. 



2.1 The N-body Simulations 

All of the simulations presented in this work were run on a 
modest Xeon cluster, using 4 nodes with dual Xeon E5520 
2.27 GHz quad-core processors per node and 24Gb shared 
memory per node. The simulations were run using the cos- 
mological structure formation software package GADGET2 
ISpringelllioOsI ). GADGET2 represents bodies by a large 
number, (in this work we use 512"^), particles. Each par- 
ticle is 'tagged' with its own unique kinematic and physical 
properties that evolve with the particle over time. GADGET2 
models the dynamics of dark matter particles using a Tree- 
PM scheme and for the purposes of this work, only dark 
matter particles were considered. 

The pre-initial particle distribution for the simulations 
used in thi s work is a g lass which has sub-Poissonian noise 
properties (|Whitelll994 ). This distribution has no preferred 
direction with forces on each particle being close to zero. If 
a glass is used as the initial condition in a standard inte- 
grator, structures do not evolve. Particle displacements are 
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imposed manually as an initial step to enable structure for- 
mation. The initial power spectrum was imposed on the par- 
ticles using the parallel initial conditions generator N-GenIC 
that was provided by Volker Springel. The initial particle 
displace ment field is form ed by using the Zel'dovich approx- 
im ation (IZerdovichlll970ll to perturb the particles, imposing 
an lEisenstehT^r^Ir l| 19981 ) matter power spectrum on the 



particles, and giving each particle an initial velocity. 

Multiple medium-resolution simulations were run with 
512'^ dark-matter particles, in a box of L = 512/1"^ Mpc 
comoving side length with periodic boundary conditions. 
The following cosmological parameters were used for a flat 
concordance ACDM mode l cons istent with the WMAP 7- 
year results (|jarosik et all I2OI0I ') : SIa = 0.73, Q.m = 0.27, 
VLb = 0.045, Us = 0.96, ag = 0.8 and h = 0.71 in units of 
100 km s"^ Mpc"^ The particle mass is 7.5 x lO^^Moand 
the softening length is 33h~^ kpc. The simulations were all 
started from a redshift of 2 = 60 and allowed to evolve to 
the present. 

The simulation data were stored at 26 output times cor- 
responding to a 128/i~^ Mpc comoving separation, between 
z = 1.5 and the present. These snapshots were chosen to fall 
within the photometric redshift error of CTz < 0.05(1 + z) 
corresponding to a displacement of ~ 147/i~^ Mpc at z — 1. 
In a 512^ particle simulation, this amounts to 100GB data 
per simulation and takes approximately 21hrs to run on the 
Xeon cluster's 32 processors. 



2.2 Shear and Convergence Map Generation 

We begin by determining the shear and convergence for a 
source plane at fixed comoving distance r^. We consider a 
distribution of sources in Section [2.41 

The effects of weak gravitational lensing on a source 
can be described by two fields, the spin-2 shear, 7, which 
describes the stretching or compression of an image, and a 
scalar convergence, k, which describes its change in angular 
size. These can be related to a lensing potential field, by 



7 = 71 + *72 = 



(1) 
(2) 



where 71 and 72 are the orthogonal components of the shear 
distortion, and d = dx + idy is a complex derivative on the 
sky. 

We want to generate shear and convergence maps along 
a lightcone through the simulation. In stead of using ray 



tracing to determine the lightcone (e.g . IWambseanss et al 
19981: iJain et al.l l200d : iTevssier et aiT l2009l : iHilbert et al 



20091 ). a line-of-sight integration was implemented using 
the Born approximatio n where one integ r ates along an 
unpe rt urbed path (e .g. 'Coorav fc Hul 120021 : IVale fc WhTt^ 
|2003| ). iFosalba et all (2008) build their convergence maps 
by adding slices from their simulation with the appropriate 
lensing weight and averaging over a pixel; 



R{e„rs) = / dr K{r,r,) S(e„r) 



(3) 



2c2 ^ 

j 



S{0.,rj) "° Ar„ (4) 

TsQj 



128h-'Mpc 




512h-'Mpc 

Figure 1. Lightcone geometry through a simulation box volume. 
The lightcone travels through the first 128/i~^ Mpc of the first 
simulation and then the next 128/i~^ Mpc of the next simula- 
tion etc. At the end of the simulation volume, the next volume 
snapshots have their centroids shifted and are randomly rotated 
to avoid repeated structures along the lightcone. 



where 6i is the position if the i* pixel on the sky and j is a 
bin in the radial direction which is at a distance of Vj and has 
a width of Arj . An overline denotes an average over a pixel 
on the sky. The expansion factor at each radial bin j is given 
by aj and the comoving radial distance of the lensing source 
plane is given by Vs . In order to make these calculations, the 
3D matter overdensity S{0, r) must be calculated by binning 
the simulation data in three dimensions. 

A limitation of this approach is memory, speed and ac- 
curacy. Here we propose, in the SUNGLASS pipeline, a new 
method for the line-of-sight integration so that no radial bin- 
ning is required to determine the convergence. The particles 
are binned in a fine angular grid while allowing them to keep 
their radial co-ordinate. 

Rewriting equation Q we find the average convergence 
in an angular pixel, with no radial binning, is given by 



K(r, 



dr K{r, 



(5) 



where AQp = AO^AOy is the pixel area and K{r,rs) is the 
scaled lensing kernel: 



K{r,rs) = 



SHpQm [r-s - r)r 
2c2 



(6) 



rsa{r) 

Hereafter we drop the overline and assume all fields are av- 
eraged over an angular pixel. A derivation of equation ((5} 
is given in Appendix|Xl In practice equation (O can be cal- 
culated by a running summation so that it is not necessary 
to re-calculate the convergence from scratch for each source 
redshift. 

The convergence maps are generated by adding the par- 
ticles that fall within the lightcone to the line-of-sight inte- 
gration. To show evolution through the lightcone, the simu- 
lation volumes are split into 128/i~^ Mpc sections. The first 
128/i~^ Mpc of the first {z — 0) snapshot is used, the second 
128h~^ Mpc of the second [z > 0) snapshot and so on until 
the end of the simulation box volume is reached at snapshot 
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Figure 2. Convergence and shear map for a simulated survey of 
100 square degrees with a single source redshift of z = 0.8. The 
colour-scale background shows the convergence while the white 
ticks show the shear signal. 



4 as shown in Figure [T] The centroid of the next simula- 
tion box is then shifted and the simulation box is rotated 
randomly to try to avoid repeated structures along th e line- 
of-sight (e.g. IWhite fc Hull2000l : IVale fc Whitll2003l ). The 
boxes are always periodic in the transverse direction. This 
continues through all of the snapshots out to a redshift of 
z — 1.5. The source redshifts have been placed at Az = 0.1 
intervals because the change in convergence between these 
redshifts is small enough that desired redshift values in be- 
tween can be accurately determined by interpolation. 

Once the convergences have been calculated at each of 
the source redshifts, the shear values can be determined on 
a flat-sky. The flat-sky shear and convergence Fourier coef- 
ficients are related by 



71 (^) 

I2(i) 



(4 



'•x \ '"y 



(7) 
(8) 



z = 0.8. There are 2048 bins in each transverse direction 
and no binning in the radial direction. The background of 
the map shows the integrated convergence along the light- 
cone up to z = 0.8 and the white ticks show the shear at this 
source redshift. The length of the ticks has been multiplied 
by an arbitrary constant to make them visible as the mag- 
nitude of the shear is at the percent level. The red patches 
show areas of the highest convergence and the shear ticks 
clearly trace these regions tangentially. These maps can be 
generated for the standard simulations at multiple source 
redshifts quite rapidly once the simulations have been run. 
The most time consuming module in this code is reading in 
the snapshots due to their reasonably large size of 100GB. 
This module can be optimised by using the fastest available 
data transfer rates on the drive where the snapshot data is 
stored. 



2.3 Shear and Convergence Power Spectra 

In order to verify the accuracy of the shear and conver- 
gence maps, the shear and convergence power spectra are 
determined for each source redshift. From equation (|4]), the 
theoretical prediction for the shear and convergence power 
spectrum for sources at redshift z is given by 



cr{z)^cr{^) 



drP 



rs{z)a'^{r] 



-,(10) 



(|Munshi et al.ll2008l ') where P{£/r;r) is the 3D matter den- 
sity power spectrum at a redshift z. 

From the simulations it is possible to determine an 
angle- averaged power spectrum from the convergence and 
shear calculated in the lightcones. When taking in to con- 
sideration the conventions used in FFTW, the discretised 
convergence power spectrum for a slice in redshift is given 
as the sum over logarithmic shells in ^-space as 



2-K 



E 



^2 Aln^' 



(11) 



where n is the total number of bins in the Fourier transform 
and A In^ represents the thickness of the shell in log £-space, 
and C'l"" is the estimated power. Similarly the shear power 
is estimated by 



where K{i) is the Fourier transform of the convergence and 



i)cr{z) 



2-K 



E 



|7i(£,z)P + |72(£,z)P 
n2 Aln£ 



(12) 



form used throughout this paper is 

FFT^M!!- 

The periodic 

nature of FFTW requires that the field is buffered with a 
small number of bins that are trimmed off after the shear 
has been calculated. To test the algorithm we also estimated 
B- modes by calculating the unphysical imaginary part of the 
convergence /3 = imag(K), from the shear, 



''X \ ''y 



— 

^x H~ ^y 



(9) 



Figure [5] is an example of a convergence and shear map 
for a field that is 100 square degrees at a source redshift of 



The Fastest Fourier 
|http://www.fftw.org/| 



Transform 



the 



West 



The B-mode power is estimated in the same way as the 
convergence. 

The modes in this power spectrum are arranged on a 
square grid, which causes discreteness errors when binned in 
annuli at small £. To correct for this, the power is scaled by 
the ratio of the measured number of modes to the expected 
number of modes, 



^max 



where g = (L/27r)^ 



(13) 



is the density of states, L is the size of 
the field in radians, £max and £min are the minimum and 
maximum wave numbers in this shell. The effect of this 
normalisation correction is about 10% at the lower wave 
numbers while the higher wavenumbers remain largely un- 
affected. The discreteness correction is not perfect which is 
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Figure 3. Simulated slices of the shear power spectra for N-body particle data at source redshifts of ^ = 0.3, 0.6, 0.8, 1.0, 1.3 and 
1.5. The smooth (red) line shows the theoretical predictions, the straight diagonal (blue) line is the predicted shot-noise at each source 
redshift. The black points are the mean power spectrum of the simulated data for the 100 realisations with errors on the mean shown 
and the (light blue) curve under these points is the simulation data with the shot-noise subtracted. The sub-shot-noise (magenta) curve 
is the estimated induced B-mode. The lower panel shows the fractional percentage difference between the simulated shear power and 
the theoretical prediction with black points representing the simulated data and light blue points representing the shot-noise subtracted 
simulation data. 
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why the same shght zig-zag of the power spectrum is evident 
in all of the source redshift planes at wavenumbers I < 100. 

We can compare our simulated shear and convergence 
power spectra with the theoretical expectation. The the- 
oretical power spectrum we use is determined using a 
code kindl y provided by Benjam i n Jo a chimi (as de mon- 



strated in I Joachimi fc Schneider! l2008l. l2009l. |201G|. an d 



extensively tested ag ainst iCosmd"l(Refregier ot al.' 120081 )'). 
This code uses the method of Smith et al. (2003) for the 
non-linear power spect rum, the matter transfer function of 
lEisenstein fc Hiil (|l998l ) and the analytical exp ression for the 
linear growth factor as given in iHeathI (| 19771 ). 

Due to the discrete number of particles in an N-body 
simulation, the measured power spectrum measured will be 
the combined real shear and convergence power plus a shot- 
noise power contribution. 



C7 



(14) 



where Cl'^ is the power estimated from the simulation. The 
shot-noise power can be derived from equation (|10|) using 
a white-noise power spectrum, Psjv(fc,r) = l/n-i{r), where 
n3(r) is the 3-D mean comoving number density of particles 
in the simulation. The shot-noise power for the shear and 
convergence is then given by 



C 



dr 



n{r)raa{r) 



(15) 



Usually, for simulated particles, n will be a constant in co- 
moving coordinates. 

Figure [3] shows the mean, normalised 2-D shear power 
spectra estimated from 100 independent simulations (black 
points and line), with the error bars showing the scatter on 
the estimated mean. The figures show the shear power for 
sources at redshifts of z = 0.3, 0.6, 0.8, 1.0, 1.3 and 1.5. 
The smooth (red) line shows the theoretical prediction for 
the ensemble-averaged shear power spectrum, while the di- 
agonal (blue) lines show the shot-noise power for each source 
redshift. The (light blue) curve between the simulated data 
and the theory curve shows the mean power spectrum with 
the expected shot-noise subtracted and the lower (magenta) 
curve shows the estimated B-mode power spectrum. 

The bottom panel of each figure shows the percent- 
age difference between the measured shear power spectrum 
and the ensemble-average theory prediction (black), while 
the lower (light blue) points show the shot-noise subtracted 
shear power spectrum. Overall the mean shear power agrees 
well, to within a few percent, with the ensemble-averaged 
theoretical model over the ^-range £ < 1000 for all source 
redshifts. The difference of a few percent is due to the fact 
that the theory 3D matter density power spectrum is a few 
percent lower than the measured data power spectrum. Cal- 
culating the highly non-linear power spectrum is currently 
not accurate to a few percent and many calculations of this 
theory curve do not agree with each other to within a few 
percent. The Joachimi theory curve was the closest fit to 
the simulations and was used for all subsequent calculations. 
At low £ the measured signal drops as we reach the size of 
the simulation box, while at high £, the estimated mean 
shear power becomes shot-noise dominated before reaching 
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^ Surveys 


Area 


n 


^median 




•2^111 ax 


100 


100 


15 


0.82 


0.05 


1.5 



Table 1. Table of mock weak lensing survey parameters used in 
this paper. 



the highest mode allowed by the resolution of the angular 
pixels beyond £ = l/0pix — lO**. 

Before reaching pixel-resolution, the measured shear 
power at high-^ agrees well with the predicted shot-noise. 
This agreement suggests that the shot-noise model works 
well in this regime , even though the i nitial particle distribu- 
tion is a glass fsee lBaueh et al.ll 19951 . for a discussion). This 
suggests an improved estimate of the mean power can be 
found by subtracting off the shot-noise contribution. How- 
ever, the shot-noise subtracted shear power does not follow 
the ensemble-averaged theoretical power estimated from the 
theory code. It is likely this is a failur e of the theoretica l 
model of lensing - on small-scales the ISmith et al.l (|2003h 
nonlinear correction formula is known to underestimate the 
matter-density power spectrum, P{k), by up to 10% at 
wavenumbers of A; < 1 and as great as 50% at fc = 10 Mpc~^ 
(Giocoli, private communication) and hence has been shown 
to underestimate the shear and c onvergence power sp ectrum 
by up to 30% on scales of ^ < 10* (jHilbert et al.ll2009l '). In the 
absence of accurate fitting formulae, simulations like those 
presented in this paper may be used to improve theoretical 
predictions. However, this needs to be explored in more de- 
tail before it is fully understood so in subsequent analysis in 
this paper we will restrict our analysis to the region of the 
measured power spectrum that agrees with the theoretical 
prediction. 

Figure |3] also shows the estimated B-mode power spec- 
tra. When galaxies trace the shear signal, we expect the 
B-mode power to pick up a shot-noise dependence. But here 
the shear signal is a pixelized field which would be con- 
tinuous in the limit of infinite pixels. Therefore we do not 
expect there to be a noise-generated B-mode. However, B- 
modes can still be generated due to leakage of power from 
the convergence field caused by the finite window function 
when we generate the shear field from equation As a 
consequence the induced B-mode has the shape of the shear 
power, but suppressed by around three orders of magnitude. 

In this section we have shown that the SUNGLASS al- 
gorithm for calculating the shear and convergence maps and 
the power spectra in redshift slices is accurate to a few per- 
cent over a wide range of scales and redshifts. Wavenumbers 
up to 1500 can be recovered for the source redshifts z ^ 1.1 
with this simulation resolution. For shot-noise subtracted 
power spectra, the recovered modes increase before the an- 
gular pixel resolution cuts off the power. 



2.4 Mock 3-D Weak Lensing Galaxy Catalogues 

Real, 3-D weak lensing data analysis is applied to a galaxy 
catalogue where galaxy angular positions and redshift are 
added to estimated shears for each galaxy. For a 2-D anal- 
ysis, individual redshifts are ignored and the theory uses 
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Figure 4. Left: The galaxy distribution, n{z), in the mock galaxy 
catalogue. The smooth (red) line shows the theoretical n{z) and 
the black histogram shows the distribution from a single simula- 
tion lightcone. The histogram shows the clustered nature of the 
lightcone. Right: The galaxy distribution in the mock galaxy cat- 
alogue with photometric redshift errors assigned to each galaxy. 
The structures visible in the true redshift lightcone have been 
smoothed out with the addition of the photo-z errors. 



only the redshift distribution. It is straightforward to gener- 
ate a simple 3-D mock weak lensing galaxy catalogue with 
the information in the lightcones we have generated from 
the simulations. Shear and convergence maps are generated 
for each lensing source redshift and then each particle in 
the simulation is assigned a shear and convergence by in- 
terpolating between adjacent planes. The error introduced 
by linearly interpolating the shear and convergence between 
source redshift planes separated by A2 = 0.1 was estimated 
by comparing with much higher redshift-sampled planes and 
found to be substantially below the theoretical prediction 
(AC/'' < 10~^) except at angular wavenumbers where shot- 
noise becomes dominant. With the interpolated shear and 
convergence assigned to each particle, we now have a fully- 
sampled 3-D mock weak lensing galaxy catalogue, which can 
be down-sampled to generate realistic weak lensing surveys. 

To down-sample the full 3-D weak lensing simu- 
lated lightcone to construct a realistic 3-D weak lensing 
galaxy catalogue, we use a galaxy redshift distribution 
(RcfrcEicr ct al. 2004) 



n{z) oc z°' exp 



(i)' 



(16) 



where zq, a and /3 set the depth, low-redshift slope and high- 
redshift cut-off for a given galaxy survey. We take o? = 2, /3 = 
2 and zo = 0.78, yielding a median redshift of Zm ~ 0.82, 
similar to the CFHTLens Survey. 

As the particles in our simulation are in comoving co- 
ordinates, we transform this redshift distribution to a prob- 
ability distribution for the particle to enter our catalogue 
given its comoving radial distance. 



pyr) Qcr I — exp 
\dzj 

where 
dr c 
" H{z)' 

and 

Hiz) = 



zo 



Ho 



[n,„(l + z)^ + Qk{1 + zY + S^a]'''' ' 



(17) 



(18) 



(19) 




1000 

Wavenumber t 

Figure 5. 2D shear power spectrum for the lightcone suite with 
the n{z) particle distribution. In the upper panel, the smooth 
(red) line is the theory prediction and the diagonal (blue) line 
is the shot noise prediction. The (black) points and line is the 
mean measured power spectrum for the suite of mock catalogues 
with the errors representing the error on the mean and the curve 
between the theory prediction and the measured simulation data 
is the shot-noise subtracted power spectrum. The diagonal (ma- 
genta) line shows the mean of the B-modes for the suite of mock 
catalogues with errors on the mean. The bottom panel shows the 
percentage difference of the data from the theory curve with errors 
on the mean (black) and the lower (light blue) points represent 
the shot-noise subtracted data. 



where Hq is the current Hubble value, fim is the current mat- 
ter density, is the current dark energy density and Q.k 
is the curvature parameter. Throughout we have assumed 
a flat, Q.K ~ 0, cosmology for our simulations. We sample 
the particle distribution so our final galaxy catalogue has 
a surface density of around 15 galaxies per square arcmin, 
with a maximum redshift cut-off at 2; = 1.5. 

The left panel of Figure |4] is an example of a redshift 
distribution taken from the full particle lightcone. The red 
line shows the theoretical distribution from eauation ll7l nor- 
malised to the number of particles selected, that the simula- 
tion particles were drawn from. The clustered nature of the 
particles in the distribution is apparent as the peaks and 
troughs around the theoretical curve can be seen. 

Our 3-D weak lensing catalogue currently assumes that 
the redshift to each galaxy is accurately known. This would 
be appropriate for a spectroscopic redshift survey, but with 
such large surveys we can expect most weak lensing cata- 
logues will contain photometric redshift estimates for each 
galaxy. To account for photometric redshift errors, we ran- 
domly sample the measured redshift from the true redshift 
using a Gaussian distribution with uncertainty 



croizg){l + Zg), 



(20) 



where Zg is the true redshift of the particle. For the purposes 
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of this work we assume a fixed tro — 0.05. The right-hand 
panel of Figure U shows what the distribution on the left 
looks like with photometric redshift errors. The structures 
are smoothed out and the distribution becomes featureless. 
The photometric redshift errors were implemented by spec- 
ifying a Gaussian error. 

Figure [5] shows the ensemble-averaged 2-D shear power 
spectrum estimated from 100 mock weak lensing surveys in 
the top panel (black dots) with errors on the mean, com- 
pared the theoretical prediction in red, and the ensemble- 
averaged B-mode power in magenta. The (blue) diagonal 
line shows the shot-noise prediction for these galaxy red- 
shift distributed lightcones. The shot-noise was determined 
by running the SUNGLASS analysis on a number of sim- 
ulation box volumes filled with randomly distributed par- 
ticles. The power spectrum of these lightcones represents 
shot-noise estimate for the simulations and is a remarkably 
straight power law. The (light blue) curve between the shot 
noise and the measured power spectrum is the shot-noise 
subtracted power spectrum. The bottom panel shows the 
fractional difi'erence between the average of the mock surveys 
and the theory curve, with the error on the mean (black) 
and the shot-noise subtracted points below (light blue). This 
shows that the mock weak lensing survey agrees with the 
theoretical expectation from wavenumbers from I = 200 to 
£ = 2000, where the disagreement with theory can be as- 
cribed to the uncertainty on the theory curve, and the rise 
of shot-noise. The shot-noise subtraction in this case is a 
few percent lower than the theoretical prediction and the 
reason for this is not well understood and is the subject of 
ongoing investigation. The analyses in this paper will use 
the measured simulation power spectrum only. The B-mode 
power appears to follow a shot-noise profile which is consis- 
tent with the efl'ect of sampling from the full particle light- 
cone. A secondary source for B-modes is source clustering, 
which appears to be sub-dominant. 

We found a dependence for the recovered shear and con- 
vergence power on the number of pixels used to estimate the 
2-D lensing power. With too many bins, there were a num- 
ber of empty pixels and this reduced the amplitude of the 
power spectrum. The amplitude of the power spectrum in- 
creased with fewer empty bins before converging at the true 
amplitude. However, by using too few bins, the number of 
I modes recovered was reduced due to pixelization effects. 
It was found that for this work, 768^ bins provided a stable 
amplitude for the power spectrum with the largest number 
of modes possible without causing this amplitude to fall. In 
this case, 0.03% of the bins are empty. If this number is in- 
creased to 5% empty, the amplitude of the power spectrum 
drops by up to 10%. This effect will also be important for ob- 
servational studies and should be considered when binning 
survey data to determine 2D lensing power spectra. 



3 PARAMETER ESTIMATION 

As described in the previous Section 100 simulations have 
been generated using the SUNGLASS pipeline. The mock 
survey parameters are given in Table [1] 

For each of these mock lensing surveys the shear and 
convergence power spectra has been estimated, and the en- 
semble average power and its scatter measured. Here we 



want to use the mock surveys to test a maximum likelihood 
cosmological parameter estimation analysis, typically used 
to extract parameters from weak lensing surveys. Here we 
try and recover the amplitude of the matter clustering, erg, 
and density parameter, f2m , from a 2-D weak lensing survey. 

In Section 12.41 we showed that our simulations could 
produce unbiased estimates of the shear power from a mock 
survey over a range of i!-modes from 200 to 2000. For pa- 
rameter estimation we need to know the conditional prob- 
ability distribution of shear power, p((77^|ct8, f^m), for the 
likelihood function, where we have fixed all other param- 
eter at their fiducial va lues. This is usually assumed to be 
Gaussian (although, see lHartlap et af]|2009l who study non- 
Gaussian likelihoods). Here we test this assumption on our 
mock catalogues. Figure [S] shows the distribution of vari- 
ations about the mean of the C^'s, AC/^, divided by the 
ensemble-averaged scatter in the power, a{C]~'). If the dis- 
tribution is Gaussian, these distributions should all lie on the 
unit-variance Gaussian. The left panel shows a histogram 
of the distribution of points for modes of ^ < 400 which 
is close to the linear region of the power spectrum. The 
middle panel shows the distribution of C]^ for modes of 
400 < £ < 1300 which represents the non-linear region of the 
power spectrum. The final panel shows the distribution for 
modes £ > 1300 which is the shot-noise dominated regime. 
The smooth (red) line in each of the panels is a normalised 
unit-Gaussian curve. In each of the panels, the histogram 
of points is peaked slightly to the left of the Gaussian peak 
which indicates a slight non-Gaussianity of the distribution 
of points. This slight non-Gaussianity may bias the Gaus- 
sian likelihood analysis but the dominant effect is currently 
the inaccurate fitt ing of the matter p ower spectrum by th e 
ISmith et al.l l|2003D formula at high k (|Giocoh et al.ll2O10D . 

The cosmological parameters of the simulations were 
estimated using Gaussian likelihood analysis where the like- 
lihood is given by 



where 



(27r)^/2(det Mu'V/^ 



exp 



-X 
2 



= }_^icr - {cr))M;,}{cj? - {cj?)), 



(21) 



(22) 



and Mill is the covariance matrix of the shear power spectra 
given by 



(23) 



The inverse covariance matrix was determined by perform- 
ing a si ngular value decom position (SVD) on the covariance 
matrix (IPress et al.lll992l l. The resulting inverse covariance 
mat rix is, however, biased due to noise in the covariance ma- 
trix. iHartlap et al] (2007) propose a correction for this bias 
by multiplying the inverse covariance matrix by a factor: 



Ns-N^ 



Ns- 



-M- 



(24) 



where Ns is the number of simulations used to determine 
the covariance matrix, Np is the number of bins in the power 
spectrum and M^J is the unbiased covariance matrix. 

The likelihood analysis relies on accurate estimation of 
the covariance matrix to show the degree of correlations. 
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I > 1300.0 




(c, - (C,)) / ■ 



(C, - (C,)) / ■ 



(C, - (C,» / cr, 



Figure 6. Histogram of the distribution of power spectra for the suite of hghtcones with the n{z) particle distribution. The left panel 
shows the distribution of the Cj~'s less than i = 400, the middle panel shows the Cj'^ distribution from 400 < I < 1300 and the right 
panel shows the distribution from ai l> 1300. 



The correlation coefficients are 

Ml,, 



(25) 



strong 



The correlation coefficient matrix is equal to 1 along the 
diagonal and the off diagonal components will show how 
correlated the I modes are, with numbers close to zero in- 
dicating low correlation and numbers close to (minus) one 
indicating high ( ant i-) correlation. 

Figure [7] shows the correlation coefficient matrix for the 
I modes being considered between 100 < I < 2500. The 
modes with a low correlation are represented in black and 
dark blues and the modes with a high correlation shown 
in yellows and reds. This shows the the bandpowers at low 
£ have very little correlation between them, as we would 
expect, since for an all-sky survey the linear power is un- 
correlated. At higher I bandpower, the modes become more 
correlated, due to cross-talk between different scales due to 
nonlinear clustering in the matter power spectrum. The vari- 
ations in this coefficient matrix indicate an error of around 
10% which is suitable for the studies in this paper. This 
error can be reduced by introducing more realisations into 
the calculations. In our analysis we shall consider modes 
up to £ = 1500, where the correlation coefficient is around 
r«' « 0.6 . 

Figure [8] shows the x^-distribution in the crg-f2m plane 
for our ensemble of simulations. The black lines represent 
the two-parameter, 1, 2 and 3(t (which should contain 
68.3%, 95.4% and 99.7% of the points assuming a bivari- 
ate Gaussian distribution), contours of parameter space for 
the cosmological parameters. However, this clearly is not 
a bivariate Gaussian distribution. The contours shown are 
representative and come from the simulation that had the 
best fit parameters that were closest to the true input pa- 
rameters (the point shown by the red polygon). The blue 
triangles represent the best fit points for each of the 100 
realisations. With this distribution, 68% of the points lie 
within the Icr contour, 93% within the 2a contour and 97% 
within the 3a. The black diamond represents the best fit for 
the combined estimate as discussed below. 

The results from this analysis give us very encouraging 
results for the parameter estimation. Figure [§] shows the 




Figure 7. Correlation coefficient matrix. This figure shows the 
correlation between the bandpower ^-modcs in the covariance ma- 
trix. The higher I bandpowers are strongly correlated (shown 
in reds), while the lower bandpowers are only weakly correlated 
(shown in blues). 



results of combining the likelihoods for all 100 realisations, 
as if we have one hundred independent 100 square degree 
surveys. Even for this test we see the maximum likelihood 
recovered parameter values lie within the 1 — cr confidence 
contour. The marginalised error on the measured parameters 
for the combined 100 surveys is Afim = 0.012 and Acrg = 
0.022, within expected errors. There is no significant bias in 
this result at the level of ~ 0.02. 



4 DISCUSSION AND CONCLUSIONS 

This work introduces SUNGLASS - Simulated UNiverses for 
Gravitational Lensing Analysis and Shear Surveys. SUN- 
GLASS is a new, rapid pipeline that generates cosmological 
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Figure 8. Gaussian likeliliood estimate. The black contours come 
from tlie simulation with the closest fit to the true cosmological 
parameters. The blue triangles show the best fit cosmological pa- 
rameters for the suite of lightcones. The true cosmological param- 
eters are shown at the red polygon and the combined best fit 
parameter is shown at the black diamond. 




0.2 



0.25 



0.3 



0.35 



Figure 9. Combined likelihood. The black lines show the 
combined Ij 2 and Scr contours. The blue triangle shows the 



best fit parameters for the combined x 
the true cosmological parameters. 



^ and the red star shows 



N-body simulations with GADGET2. It computes weak lens- 
ing effects along a lightcone using line-of-sight integrations 
with no radial binning and the Born approximation to deter- 
mine the convergence and shear at multiple source redshifts. 
This information is interpolated back on to the particles in 
the lightcone to generate mock shear catalogues in 3D for 
testing weak lensing observational analysis techniques. 

In this work, SUNGLASS was used to generate 100 sim- 
ulations with 512'^ particles, a box length of bl2h~^ Mpc 
and a WMAP7 concordance cosmology. The corresponding 
mock shear catalogues were 100 sq degrees with a source red- 
shift distribution with median Zm = 0.82 and 15 galaxies per 
square arcminute. The parameters are easily changed within 
the SUNGLASS pipeline so that the mock shear catalogues 
matches the survey of interest. 

To show the reliability of the lightcones generated with 
SUNGLASS, E- and B-mode power spectra were shown at 
multiple source redshifts. The results show that at low red- 
shifts, the signal becomes dominated by shot-noise at rea- 
sonably low I. With increasing source redshift, the power 
spectrum recovers the theoretical prediction over a wider 
range of modes, I < 2500. 

Given that the measured power spectrum of the simu- 
lations appears to follow the predicted shot noise at higher 
modes, the shot noise was subtracted from the power spectra 
to increase the recovered range. The theoretical prediction 
is expected to under predict the power spectrum around the 
turn over and consequently, the simulations could be recov- 
ering the power spectrum up to around £ — 5 x 10* at the 
highest redshift planes. 

The multiple source redshift plane shear and conver- 
gence was interpolated onto the particles in the lightcone to 
generate a mock shear catalogue. A redshift sampling was 
also imposed on the lightcone to mimic an observed shear 
catalogue. Binning this distribution too finely resulted in 
empty bins which had the effect of suppressing the power 
spectrum. This has implications for observations where the 
number of objects per square arcminute should be taken 
into account, as well as the density of the binning, when 
determining the accuracy of the power spectrum. 

The mock shear catalogues were used to determine a 
covariance matrix which is essential for both parameter es- 
timation and data analysis. A strength of SUNGLASS is the 
ability to rapidly produce Monte Carlo realisations of these 
catalogues, ensuring independent mock data sets for the gen- 
eration of the covariance matrices. 

The mock catalogues were also used to perform a sim- 
ple parameter estimation using Gaussian likelihood analysis. 
The distribution of power spectra were shown to be reason- 
ably Gaussian and the resulting parameter estimation con- 
tours for a single realisation showed a good agreement with 
the input parameters within the 2-parameter 1,2 and 3a er- 
ror contours. 

The combined likelihood from the 100 simulations 
shows narrow likelihood contours and accurate parameter 
recovery within the expected errors, with no evidence of sig- 
nificant bias at the level of ~ 0.02. 

Current and future telescope surveys promise to pro- 
vide an enormous amount of data for weak lensing analysis. 
Weak lensing is still a young field and analysis techniques are 
still being developed. It is essential that the strengths and 
weaknesses of these techniques are fully understood before 
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using them on real data with unknown parameters. Using 
the simulations, lightcones and mock shear catalogues pro- 
vided by the SUNGLASS pipeline, and demonstrated in this 
paper, is an excellent way to test these observational weak 
lensing analysis techniques. The outputs of this pipeline 
have been rigorously tested and are well understood, mak- 
ing them ideal for generating covariance matrices that are 
critical to many observational analysis techniques. 
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APPENDIX A: DERIVATION OF THE 
LINE-OF-SIGHT CONVERGENCE WITH NO 
RADIAL BINNING 

This appendix shows how the line-of-sight convergence 
shown in equation [S] was derived. 

Start with the general equation for the convergence. 



AC = / dr K{r,rs) 5{r), 



(Al) 



where Vs is the lensing source redshift, 5{r) is the fractional 
matter overdensity and K{r,rs) is the kernel 



(r^3iW._ 
rsa[r) Zc-' 

The overdensity 5{r) is given by 
n{r) 



(A2) 



(A3) 



where n(r) is the average density at the comoving radial 
distance r and is constant in comoving co-ordinates. 

The particle number density, n(r), is given by a sum of 
3D delta functions 



^^^5^-(0-0O,(A4) 
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where part are the particles in the pixel with 55 Ts. Substi- 
tuting this sum of delta-functions into equation (|A1|) yields 
the average convergence per pixel on the sky, p, with no 
radial binning; 




where Ar2p = AO^Ae, 
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